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Abstract 

Most aquatic vertebrates swim by lateral napping of their bodies and caudal fins. 

^U| While much effort has been devoted to understanding the flapping kinematics and its 

influence on the swimming efficiency, little is known about the stability (or lack of) 
of periodic swimming. In this paper, we examine the stability of periodic locomotion 
due to sideways flapping in unbounded potential flow. It is believed that stability 
limits maneuverability and body designs/flapping motions that are adapted for stable 
swimming are not suitable for high maneuverability and vice versa. Here, we consider 
a simplified model where the swimmer is a planar elliptic body undergoing prescribed 
periodic heaving and pitching. We show that periodic locomotion can be achieved 
due to the resulting hydrodynamic forces, and its value depends on several parameters 
O including the aspect ratio of the body, the amplitudes and phases of the prescribed 

flapping. We obtain closed-form solutions for the locomotion and efficiency for small 
flapping amplitudes, and numerical results for finite flapping amplitudes. We then 
study the stability of the (finite amplitude flapping) periodic locomotion using Floquet 
theory. We find that stability depends nonlinearly on all parameters. Interesting 

CNJ trends of switching between stable and unstable motions emerge and evolve as we 

continuously vary the parameter values. This suggests that, when it comes to live 
organisms, maneuverability and stability need not be thought of as disjoint properties, 
rather the organism may manipulate its motion in favor of one or the other depending 

**) on the task at hand. 

o 
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A large proportion of fish species are characterized by elongated bodies that swim forward by 
flapping sideways. These sideways oscillations produce periodic propulsive forces that cause 
the fish to swim along time-periodic trajectories, [I]. The kinematics of the flapping motion 
and the resulting swimming performance, as well as their relationship to the swimmer's 
morphology, have been the subject of numerous studies, see, for example, [21 E]. However, 
little attention has been given to the stability of underwater locomotion. The importance 
of motion stability and its mutual influence on body morphology and behavior is noted in 
the work of Weihs, see [I] and references therein. In [3], Weihs used clever arguments and 
simplifying approximations founded on a deep understanding of the equations governing 
underwater locomotion to obtain "educated estimates" of the stability of swimming fish 
without ever solving the complicated set of equations. 



The main objective of this paper is to study the stability of periodic locomotion result- 
ing from periodic sideways flapping. The locomotion is said to be unstable if a perturbation 
in the conditions surrounding the swimmer's body result in forces and moments that tend 
to increase the perturbation, and it is stable if these emerging forces tend to reduce such 
perturbations or keep them bounded so that the fish returns to or stays near its original 
periodic swimming. 

Stability may be achieved actively or passively. Active stabilization requires neurolog- 
ical control that activate musculo-skeletal components to compensate for external perturba- 
tions acting against stability. On the other hand, passive stability of the locomotion gaits 
requires no additional energy input by the fish. In this sense, one can argue that stability re- 
duces the energetic cost of locomotion. Therefore, from an evolutionary perspective, it seems 
reasonable to conjecture that stability would have a positive selection value in behaviors such 
as migration over prolonged distances and time. However, stability limits maneuverability, 
and body designs/flapping motions that are adapted for stable swimming are not suitable 
for high maneuverability and vice versa, [U [5] . 

In this work, we study stability of periodic swimming using a simple model consisting 
of a planar elliptic body undergoing prescribed flapping motion in unbounded potential flow. 
By flapping motion, we mean periodic heaving and pitching of the body as shown in Figure [TJ 
We formulate the equations of motion governing the resulting locomotion and examine its 
efficiency We then investigate the stability of this motion using Floquet theory (see [6]). We 
find that stability depends in a non-trivial way on the body geometry (aspect ratio of the 
ellipse) as well as on the amplitudes and phases of the flapping motion. Most remarkable is 
the ability of the system to transition from stability to instability and back to stability as 
we vary some of these parameters. 

This model is reminiscent of the three-link swimmer used by Kanso et al. to examine 
periodic locomotion in potential flow, see [7]. The three-link swimmer undergoes periodic 
shape deformations that result in coupled heaving, pitching and locomotion. Here, we ignore 
body deformations for the sake of simplicity and prescribe the heaving and pitching motion 
directly. Note that the three-link swimmer was also used by Jing & Kanso to study the effect 
of body elasticity on the stability of the coast motion offish (motion at constant speed). They 
found that elasticity of the body may lead to passive stabilization of the (otherwise unstable) 
coast motion, see [HI IS]. The present model consisting of a single elliptic body is mostly 
similar to the system studied by Spagnolie et al. both experimentally and numerically, [10J. 
In the latter, an elliptic body undergoes passive pitching (via a torsional spring) subject to 
prescribed periodic heaving in viscous fluid, whereas in our model both the pitching and 
heaving motions are prescribed and the fluid medium is inviscid. Despite these differences, 
the two models exhibit qualitatively similar behavior as discussed in Section [6] 

The paper is organized as follows. In Section [2j we formulate the equations of motion 
governing the locomotion of a periodically flapping body in unbounded potential flow. We 
analyze the body's locomotion and efficiency when subject to small amplitude flapping mo- 
tion in Section [3j and consider the more general case of finite amplitude flapping in Section |4j 
In Section |5j we assess the stability of the periodic locomotion using Floquet theory. The 
main findings and their relevance are discussed in Section [6j 
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Figure 1: Model of flapping fish: An ellipse with semi-axes a and b is submerged in unbounded potential 
fluid. Motion is observed in inertial frame with position of mass center given by (x, y) and orientation by 9. 
Fish flaps in y- and 0-directions, and propels in x-direction. 

2 Problem Formulation 

Consider a planar elliptic body with semi-major axis a and semi-minor axis b, submerged in 
unbounded potential flow that is at rest at infinity. The elliptic body is neutrally buoyant, 
that is to say, the body and fluid densities are equal to p. The body's mass is given by 
m b = pnab, and its moment of inertia about the center of mass C is J b = m b (a 2 + b 2 )/4. Let 
(x, y) denote the position of the mass center with respect to a fixed inertial frame and let 
denote the orientation angle of the ellipse measured from the positive x-direction to the 
ellipse's major axis, see Figure II} The linear and angular velocities are given by (x, y) and 
0, respectively, where the dot () correspond to derivative with respect to time t. 

In order to emulate the flapping motion of a swimming body, we assume that y and 9 
vary periodically in time due to some periodic flapping force F flap and flapping moment r flap 
generated by the swimming body. Note that in the case of a body swimming by deforming 
itself, y and 9 are a result of the body deformation. Here, we do not account for the body 
deformation but rather prescribe y(t) and 9{t) as periodic functions of time. Namely, we set 

y{t) = A y sm(ut + (j)y), 9{t) = A e sin(wt + <f> 9 ). (1) 

and solve for the resulting locomotion in the x-direction. 

The equations governing the motion of the flapping body are basically Kirchhoff's 
equations expressed in inertial frame and subject to forcing _F flap and r flap in the y- and 
^-directions, that is, 

m b x = F x , (2) 

and 

m b y = F y + F flap , J b 9 = r + r flap , (3) 

where F x , F y and r are the hydrodynamic forces and moment acting on the body. For motions 
in potential flow, F x , F y and r can be obtained using a classic procedure, see Appendix A, 

1 1 

F x = - [—(mi + 1712) + \JTi2 — TTi-i ) cos 29] x H — (rri2 — mi)y sin 20 — (777,2 — m>i) {x sin 29 — y cos 29)9, 

1 1 

F y = - [—(7721 + 1712) — (777.2 — 777,1) cos 29] y -\ — (7772 — mi)xsm29 + (777,2 — mi) (i cos 2(9 + y sin 29)9, 



r = —JO + -(7772 — 777i) (x 2 sin26 ) — y 2 sin.20 — 2xycos20) 



(4) 



Here m x = pnb 2 , m 2 = pna 2 are, respectively, the added mass of the elliptic body along its 
major and minor directions and J = pii(a 2 — b 2 ) 2 /8 is the added moment of inertia [TT] . 
Substituting Q into ^ and ([3]), one can use §2§ to solve for x(t), and (|3| to compute the 
forcing F p and r p needed in order for the body to achieve the prescribed flapping motion 
in©. 

The total angular momentum h of the body-fluid system is given by h = (Jb + J)9 
whereas the total linear momentum can be written as the following 

p x \ _ fi\ /cos# — sin^X (m\ \ / cos# sin^X (x\ , , 

PyJ \yJ \sin6 1 cos 9 / \ m 2 J \— sin9 cos 9 J \yj 

The momentum p x is conserved since there is no external forcing applied in the x direction. 
Therefore, one has p x (t) = p x (t = 0), which yields 

*(*) = of ~ mi ~ 2 Tl rr^ [(y sin 29) t=0 - (y sin 29)] . (6) 

2 [nib + Wi cos z o + m 2 sin 1 

That is to say, equation (|2| admits an integral of motion whose value is given by the above 
equation. The distance traveled by the body's center of mass in one period of flapping, 
T = 2ti/uj, is given by 



d 



r 

xdt 
o 



\x(T)-x(0)\. (7) 



The total kinetic energy E of the body-fluid system is given by 

E= l -(xp x + yp y )+ l -(J b + J)9 2 . (8) 

By the work-energy theorem, the time derivative of the kinetic energy is equal to the total 
power input by the flapping force F flap and moment r flap . Thus, the work done by flapping 
is equivalent to the kinetic energy E. To this end, the average work done in one period is 
given by 

E = ±f E(t)dt. (9) 

We define the cost of locomotion e as the average work divided by the average distance over 
one period, namely 

e = f. (10) 

Hence, smaller e means less energy expenditure for a fixed distance traveled. Therefore, it is 
convenient to denote the efficiency of the system 7] as the inverse of the cost of locomotion, 
that is r] = 1/e = d/E. 

Before we proceed to examining the locomotion and efficiency of such swimmer, it is 
convenient to non-dimensionalize the system by scaling time with T, length with yah and 
mass with nib — pnab. The variables are subsequently written in dimensionless form. The 
important parameters for this system are: aspect ratio 7 = a/b, amplitudes A y and Ag, and 
phases (fi y and (fig. 



3 Small Flapping Amplitudes 

Consider the case with small flapping amplitudes A y and Ag, that is to say, let A y = e y <C 1 
and Ag = eg <^ 1 where both e y and eg are of the same order of magnitude. Hence, one 
has y,y,y ~ 0(e y ) and 9,8,9 ~ 0(eg), but u, cfi y and (fig are not necessarily small. Take the 
approximation cos 9 « 1 and sin 9^9 and substitute into (pi) to obtain 



^ ~ o (7 - 1) we^eg [sin(2wt + y + e ) - sin((fi y + e )] . 



Ill 



Clearly, the velocity in the x direction depends on the aspect ratio 7 and (fi y + (fig. This 
suggests that as long as (fi y + (fig = 2nir + constant, x is the same function of time. Its 
magnitude is of order ~ 0(e y eg). In other words, for small amplitude flapping, the motion in 
x is small compared to the flapping motion in y and 9. Approximating expressions of F flap 
and r flap are obtained by substituting (fTl) and (11) into pi), which yields 



F flap w (m b + m 2 )y, r flap « ( J h + J)9. (12 

For small amplitude flapping, we can express the cost of locomotion in closed form 

vr [7(7 + 1)4 + 2 (7 2 + 1)4] 



27(7 - l)e y e e \ sin((f)y + (f)g)\ ' 
Hence, to minimize e (or, equivalently, to maximize efficiency if), one needs 

1' 



(13) 



7 — y 00, 



and 



+ 



n + 



TV. 



n 



0, ±1, ±2, 



(14) 



For large amplitudes A y and Ag, one does not have a closed form expression for e, hence the 
efficiency needs to be analyzed numerically, as done in the next section. 

4 Locomotion and Efficiency 

In this section, we examine the swimming trajectories and their dependence on the following 
parameters: aspect ratio 7, amplitudes A y and Ag, and phases (fi y and (fig. The swimming 
motion is given by ([I]) and ^, where the latter is integrated numerically to get x(t). 
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Figure 2: (color online) Left: Prescribed napping in (y,0) plane. Initial points are marked by o. Right: 
Trajectories of mass center in (x, y) plane with snapshots of body in motion overlaid. Simulations are for 
A y = 1, Aq = f , 4> y = — \,4>g = and various aspect ratios: (a) 7 = 1.01, (b) 7 = 4, (c) 7 = 8, (d) 7 = 1000. 
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Figure 3: (color online) Left: Prescribed napping in (y,0) plane. Initial points are marked by o. Right: 
Trajectories of mass center in (x, y) plane. Simulations are for Ag = ? , 7 = 4, <fi y = — ^,<fig = and various 
heaving amplitudes: (a) A y = 0.01, (b) A y = 0.5, (c) A y = 1, (d) A y = 2. 
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Figure 4: (color online) Left: Prescribed flapping in (y,0) plane. Initial points are marked by o. Right: 
Trajectories of mass center in (x, y) plane. Simulations are for A y = 1,7 = 4, y = — f , 06> = and various 
pitching amplitudes: (a) Ag = 7r/8, (b) Ag = n/4, (c) Ag = tt/2, (d) Ag = 3n/4. 
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Figure 5: (color online) Left: Prescribed flapping in (y,0) plane. Initial points are marked by o. Middle 
and right: Trajectories of mass center in (x, y) plane. Simulations are for A y — 1, Ag = j, 7 = 4 and various 



combinations of phases: (a) (<j}y,(j>e) = (— f,0), (b) 



(" 



3tt 



"y,> 



= (-f.-f). (C) (^,^) = (f ,0), (d) (0 V ,M = 



), (e) (0„,M = (-f,0), (f) (^<t>g) = (0,0), (g) (0,> e ) = (-f,-f), (h) (^,</, e ) = (f,-2f ). 



Consider the case where A, = 1, Ai 



4' Vy 



— f > 06* = and consider various aspect 



ratios 7 = 1.01,4,8, 1000, as shown in Figure |2j Note that as we vary the aspect ratio, the 
total area of the elliptic body is maintained constant, which is guaranteed by the way we 
non-dimensionlize length using yob. As expected, the net locomotion is almost zero when 
the elliptic body is close to a circular shape (7 = 1.01) and it reaches a maximum as the 
elliptic body approaches a flat plate (7 = 1000). 

In Figure [3j 7 is set to 4 and A y is varied. One can see that the net locomotion d 
depends linearly on A y , which is also evident from ([6]). In Figure |4l different cases of Ag are 
shown. Here, the net locomotion depends nonlinearly on Ag. Interestingly, the trajectories 



that correspond to Aq = it/A and A$ = n/8 almost coincide, whereas for Aq = 37r/4 the 
locomotion is in the negative x direction. 

Motions for various phases <p y and <pg are shown in Figures [5] Notice the shape and 
orientation of the closed path in the (y, 9) parameter space depend on the difference in phase 
4> y — (pe. This can be readily verified by eliminating t from (fl| and expressing the closed 
path in the (y, 9) plane as 

i$- 2 ii^-^ + {iS= si ^'- M - (15) 

In other words, as <fi y — <pg varies, the closed path in the (y, 9) plane is elliptic, except for 
<Py — <Pe — niT {n — 0, ±1, ±2, . . .) in which case it is a segment of the straight line given by 
9 = (—l) n (Ag/A y )y. Further, notice that the the dependence of resulting locomotion x(t) 
in d6]) on the flapping motion x(y(t),9(t)) possesses the following symmetries 



x(-y,-0) =x(y,9), x(-y,9) = -x(y,6), x(y,6) = -x(y,6). (16) 

whereas the flapping motion in (II]) has the following symmetries 



y(t; <p y + tt) = -y{t; </> y ), y(-t; <f> y ) = -y(t; -<j> y ), 
9{t- <f) d + 7i) = -9{t; fo), 9{-t; fo) = -9{t; -</» e ). 



(17) 



Based on these symmetries, one can immediately conclude that, when all other parameters 
are held fix, motions that correspond to (<f> y , 4>e) an d (—(j> y , —<fie) are mirror images of each 



other: their distances and energies are the same, which can be seen from (16). Also, for 
combinations (<p y ,<pe) with cf) y — <ftg = 2mc + constant, they result in the same path in the 
(y, 9) space. When tracing the same path in the (y, 9) plane (but starting at different initial 
points), the resulting trajectories in the (x,y) plane are very similar (with different initial 
positions). One might be tempted to jump to the conclusion that tracing a straight line 
in the (y, 9) plane corresponds to zero net locomotion in the (x,y) plane. This is true for 
4>y = 4>e = and <j) y = (pg = — tt/2, but it is not true in general as evident from the example 
of 4> y = 7r/4, 4>e = — 37r/4 in Figure |5F/i). This is unlike the well known result in geometric 
mechanics that zero enclosed area in shape space results in zero net locomotion, see [7]. 
Indeed, the locomotion here is not a result of a geometric phase but rather of a dynamic 
phase. 

While the numerical results in Figures [2] to [5] provide an understanding of the average 
distance d travelled over one period T of flapping, the average work E is not reflected and 
remains to be computed. More specifically, we are interested in the cost of locomotion e = 
E/d, which we compute numerically. Ideally, one would like to find optimal parameter values 
that minimize e (maximize efficiency 77) and/or maximize d (see, for example [31 H2J). It is 
computationally expensive to do an exhaustive search to minimize e over the five dimensional 
parameter space. Instead, similar to the locomotion study, we focus on varying one parameter 
at a time. In Figure [6J we set (f> y = — |, (f)g = while 7, A y and Aq are varied in (a), (b) and 
(c), respectively. Solid lines are nonlinear numerical solutions, and dashed lines are obtained 
by substituting the parameters into (13). Figure [61(a) shows that, for A y = l,A e = n/A 



there exist a optimal value of 7 ~ 2.9, whereas the small amplitude approximation in (13) 
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Figure 6: Cost of locomotion e as a function of: (a) aspect ratio 7, (b) heaving amplitude A y , (c) pitching 



amplitude Ag. The base parameter values are set to 



-7r/2, <^ e = 0, 7 = 4, Ag = tt/4, A, = 1. Solid 



lines are nonlinear numerical solutions, while dashed lines are based on small amplitude approximation given 



in ( 13 ) 



(shown in dashed line in the figure) predicts that e is a decreasing function of 7. Similarly, 
Figures Mb) and (c) show optimal values of A y ps 1.3 and Ag ps ||. One can see that the 
small amplitude results qualitatively follows the nonlinear behavior of e for the parameter A y . 
This is not surprising given that the work E depends quadratically on A y , the displacement 



d depends linearly on A y and that the small amplitude approximation in (13) preserves 
the form of dependence on this parameter. However, in the case of varying Ag, the small 
amplitude results can approximate the nonlinear efficiency up to Ag w ||, but start to 
deviate from the nonlinear results thereafter. 

In Figure [7] to |9l we examine the dependence of e on the phases by discretizing ((j) y , <pg) 
on a 201 x 201 mesh in the range of [— n , n] x [— n , n] (which is representative of the entire 
plane due to periodicity). Contours of e are depicted for various 7, A y and Ag as a function 
of (4> y ,4>g). Indeed, it is sufficient to only show the dependence of e on one quarter of 
the shown region, say, [0 , 7r] x [0 , n], because of the reflection symmetry about the origin 
e (</>„, 4>o) = e(—<j>„, —4>e) and the periodicity of n in efficiency e(<fi y + n, 4>g + ir) = e( ' 



'y, 



'y: 



both of which can be easily obtained from (16 17). 



J yi 



Note that the parameters that minimize e, thus maximize efficiency 77, are approxi- 
mately 3 < 7 < 4, A y ?z 1.3, Ag ?z 37r/16 and (4> y ,(j)g) ~ ((m + ^)rr,mr), where m,n — 
0, ±1,±2, ... One example of the optimal phases is 4> y = ^,<f>o — 0, whose corresponding 
locomotion and snapshots of dynamics are shown in Figure [2} For this optimal motion, 
the pitching angle is zero when the heaving motion is maximum (90° out of phase), which 
qualitatively agrees the results in [12] . The optimal aspect ratio 3 < 7 < 4 agrees with the 
optimal shape aspect ratio obtained in the optimization study in [3], and is representative 
of the aspect ratio of various Carangiform swimmers such as bass (7 = 3.8) in [13], tuna 
(7 = 3.5) in [13] and saithe (7 = 4.1) in [15]. The optimal heave to cord ratio A y /a ~ 0.75 
(where a = */7 ~ \/3) and maximum angle Ag m 37r/16 = 16.875° both agree with the 
optimal motions for the rigid flapping body (0.75 < A y /a < l,Ag ps 16°) given in 



(a) 7 = 1.01 



(b) 7 = 2 



(c) 7 = 4 



(d) 7 = 8 



(e) 7 = 16 



Figure 7: Contour plots of cost of locomotion e for the cases A v — 1,A$ = ir/4 and various aspect ratio 
7. Each plot is evaluated on a 201 x 201 mesh in [—it, tt] x [— tt , tt] in {(j) y ,(j)g) plane. Lower value of e 
corresponds to higher efficiency. 




(a) A y = 0.01 



(b) A v = 0.5 



(c) A y = 1 



(d) Ay = 2 



(e) A y = 6 



Figure 8: Contour plots of cost of locomotion e for the cases Ag — 7r/4, 7 = 4 and various heaving 
amplitude A y . Each plot is evaluated on a 201 x 201 mesh in [—tt, tt] x [—tt, tt] in (<f) y ,(f)0) plane. Lower 
value of e corresponds to higher efficiency. 




(a) A e = 0.01 



(b) Ag = TT/ 



(c) Ag = tt/4 (d) Ag = 3^/8 (e) Ag = tt/2 



Figure 9: Contour plots of cost of locomotion e for the cases A y = 1, 7 = 4 and various pitching amplitude 
Ag. Each plot is evaluated on a 201 x 201 mesh in [—TV, tt] x [—tt, tt] in (4> y ,(f>g) plane. Lower value of e 
corresponds to higher efficiency. 



5 Stability of Periodic Locomotion 

We study stability of periodic motion subject to arbitrary perturbations in the surrounding 
fluid. We begin by introducing q = [x, y, 9, 9} T , and rewriting ^ and ^ as follows 



flap 



q = f (q) + F 

where detailed expressions for M, f and F flap are listed in Appendix B. In Sections |2]-[4j we 
prescribed the flapping motion y(t) and 9(t) according to ([l]) and used ^ to solve for x(t) 
and ([3]) to solve for F flap and r flap necessary to produce the flapping motion. The resulting 
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Figure 10: Left: stability of the case A y = 1, Ag = tt/4 and 7 = 4, with phases (<j> y ,4>e) varied on a 
201 x 201 mesh in [— tt , 7r] x [— 7r , 7r]. Stable cases correspond to shaded areas, unstable cases are white 
areas. Middle and right: (top) (<j> y , <pg) = (—tt/2, 0) is stable. Solid lines correspond to unperturbed periodic 
solutions, dashed lines correspond to perturbed solutions; (bottom) (</> y ,(j>g) = (— tt/4, — tt/4) is unstable. 



motion x(t), y(t), and 9(f) as well as the forcing F flap and r flap are periodic with period T. 
We let q p denote the q corresponding to such periodic motion. We study the stability of q p 
by introducing a small perturbation <5q such that q = q p + 5q while keeping _F flap and r flap 
the same as that producing the periodic solution. In other words, we account for arbitrary 
perturbations in the fluid environment while keeping the same flapping forces to check if 
such perturbations destabilize the periodic trajectory. 



Linearizing equation (18) about the periodic trajectory q p (£), one gets 

<5q = J(t)<Jq, 
where the Jacobian J(t) is a 4 x 4 matrix, and it has period T, that is, J(0) = J(T), 



S(t) 



dg 
<9q 



where 



(qp,Fflap) 



_1 (f + F 



flap^ 



(19) 



(20) 



See Appendix B for details. Let $(£) denote the fundamental solution matrix of (19). The 
eigenvalues Aj of the time-independent matrix, 



B = $(0) _1 $(T). 



(21) 



are referred to as the characteristic multipliers. Their locations in the complex plane indicate 
the stability of the periodic solution q p : if any of the Aj's lies outside the unit circle, then q p 
is unstable; if all A, are on the unit circle, then q p is regarded to as marginally stable. For 
a non dissipative system like our model, these two are the only possible scenarios, namely, 
unstable or marginally stable (simply referred to as "stable" hereafter). One eigenvalue Ai 
is always 1, reflecting the fact that q p is periodic. The other eigenvalues may be complex. 
If, say, A2 and A3 are complex, they come in conjugate, i.e. A2 = A3, the fourth eigenvalue 
A4 must be real. 

We determine stability by computing the characteristic multipliers of the 4-dimensional 
system described in (19). For the case A y = l,A 9 = |,7 = 4, the stability results are 
plotted in Figure 10 as a function of the phases ((f> y ,(f>g), again evaluated on a 201 x 201 
mesh in [— tt , tt] x [—tt, tt]. The stable regions are shaded areas, and the unstable regions 
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Figure 11: (Color online) Characteristic multipliers of the cases A y = l,Ag = 7r/4,7 = 4, <j) y = 0, and 4>e 
varied from — 7r to n: (a) Real and Imaginary parts in complex plane. Two characteristic multipliers always 
locate at (1,0) are represented by *. The two complex conjugates are represented by O- 3 ) R- e& l an d ( c ) 
Imaginary parts of characteristic multipliers. Motion becomes unstable when the two complex conjugates 
collide at (1, 0) and split on real axis. 



are white areas. Notice the reflection symmetry about (0, 0) and periodicity of 7r in both cf) y 
and <j)g that we observed in the efficiency analysis is again seen in the stability plot. Two 
examples with different stability characteristics are shown. The solid lines are unperturbed 
periodic solutions q p , and dashed lines are solutions with random initial perturbations with 
magnitude |<5q(£ = 0)| ~ O(10~ 3 ). Clearly, the trajectory corresponding to parameters in 
the stable region remain close to the periodic trajectory for the integration time whereas 
that corresponding to parameters in the unstable region does not. 

In Figure ITT1 we examine the behavior of the real and imaginary parts of the charac- 



for A y = l,Ag = vr/4, 7 = 4 and 



0. In other 



J yi 



space along the dashed line 



teristic multipliers as a function of ry n,i ^ y 

words, we explore the behavior of Aj's as we cross the 

corresponding to <p y = in the left plot of Figure 10 One can see that two characteristic 

multipliers are always located at (1,0) (represented by it). The other two are represented 

by O- The dynamics changes from stable to unstable when the two conjugates collide at 

(1, 0) and split onto real axis. For the considered parameters, when <pg varies from — it to n, 

stability changes from unstable to stable and stable to unstable four times in total. 

We now show how stability changes with 7, A y and Ag by varying them one at a time 
while keeping the others constant. 



l,A e 



Figure 12 shows stability regions for A y — .... . ,,., - 
and various aspect ratio 7. For bodies closer to circular shape (7 = 1.01), the motion is 
stable for all ((f> y , <pg) (but not useful since net displacement is almost zero). When 7 > 1.43, 
unstable regions start to appear as strips. As 7 increases, unstable regions grow while stable 
regions shrink. The total area of stable regions becomes minimum when 7^2, and the 
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stable 



(a) 7 = 1.01 






(b) 7 = 1.45 



(c) 7 = 1.75 



(d) 7 - 2.5 




(e) 7 = 4 



(f)7 = 



(g) 7 = 16 



(h) 7 = 32 



Figure 12: Stability regions for the cases A y = l,Ag = 7r/4 and various aspect ratio 7. Each plot is 
evaluated on a 201 x 201 mesh in [— 7r, it] x [—w, it] in (4> v ,(J)q) plane. Shaded areas correspond to stable 
cases, white areas correspond to unstable cases. 



stable areas remain around (4> y ,(f)g) = (( m + l/2)7r, mr). Interestingly, as 7 continues to 
increase, new stable regions start to emerge and grow from the previous unstable areas 
around (cf) y ,4>g) = (mr, (m+ l/2)7r). Then, at these spots, unstable regions emerge and grow 
from the newly formed stable regions, and so on and so forth. The boundaries between 
stable and unstable regions around (rnr, (m + l/2)7r) become blurry as 7 becomes larger, and 
((m + l/2)7r, rnr) remain stable. Overall, the total area of stable regions recovers. This trend 
is reminiscent of the phenomenon observed in Spagnolie et al. [10], in which the authors 
noticed the motion of an elliptic body subject to prescribed heaving motion and passive 
pitching motion goes through states from "coherence to incoherence, and back again" as the 
aspect ratio changes. Note that the latter studies are in viscous fluid whereas the analysis 
here is for an inviscid fluid model. Still, this simplified model is able to capture, at least, 
qualitatively, the behavior observed in [TO] . 

Figure 13 shows stability regions for Ag = 7r/4,7 = 4 and various A y . For small A y 
such as 0.01, the body is mostly rotating, and the motion is stable for all (4>„, <t>g) but with no 



'yi 



net locomotion. As A y increases, unstable regions start to form around (nir, (m + l/2)7r), as 



can be seen in Figure |13[c). As A y continues to increase, unstable regions grow while stable 
regions shrink. Then, layers of stable/unstable regions start to form around (nir, (m + l/2)7r). 
Unlike previous trend when 7 is varied, areas around ((m + 1/2)tt, wk) do not remain stable. 
Instead, the stable/unstable layers also form there. Overall, the total area of stable regions 
decreases as A y becomes larger. The area and shape of stable regions depend nonlinearly on 
A y , which is an interesting result since the trajectory of the mass center depends linearly on 

Ay. 

Finally, we vary Ag while keeping A y = 1 and 7 = 4 in Figure 14 When Ag is small, 
the whole plane in (<p y , (pg) is stable but again does not result in net locomotion. As Ag 
increases, unstable regions start to emerge and grow, while stable regions shrink and remain 
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4>e 



(a) A y = 0.01 (b) A y = 0.6 (c) A y = 0.625 (d) A y = 0.75 



(e) A y = 1 



lOt 

(f) A a = 1.25 



),UJi 
me 







<*„ 




(g) A y = 1.5 (h) A y = 1.75 



(i) Ay = 2 



(J) Ay = 6 



Figure 13: Stability regions for the cases Ag — ir/4, 7 = 4 and various heaving amplitude A y . Each plot 
is evaluated on a 201 x 201 mesh in [— -k , n] x [— n , 7r] in {(j) yi <fig) plane. Shaded areas correspond to stable 
cases, white areas correspond to unstable cases. 



stable 



ft A 




(a) A e = 0.01 



(b) Ag = tt/8 (c) Ag = 371-/16 (d) Ag = 7r/4 




(c) Ag = 5tt/16 (f) Ag = 3tt/8 (g) Ag = 7tt/16 (h) Ag = ir/2 

Figure 14: Stability regions for the cases A y = 1, 7 = 4 and various pitching amplitude Ag. Each plot is 
evaluated on a 201 x 201 mesh in [— ir , n] x [— w , ir] in (4> y ,(j)g) plane. Shaded areas correspond to stable 
cases, white areas correspond to unstable cases. 

around (mr, (m + 1/2)tt). However, as Ag continue to increase, the stable regions gradually 
become strip-like shaped roughly along (fig — <f) y = {m + |)7r. And unstable regions start to 
emerge within the stable strips. 

This trend of switching from stability to instability and back to stability that is observed 
when varying each of the parameters is extremely interesting. It suggests that such swimmers 
can change their stability character by changing their flapping motion, and thus can easily 
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switch from stable periodic swimming to an unstable motion (more maneuverable) when 
they feel the need to, such as when evading a predator. Based on this, one can conjecture 
that when it comes to live organisms, maneuverability and stability need not be thought of 
as disjoint properties, rather the organism may manipulate its motion in favor of one or the 
other depending on the task at hand. Experimental results that show that live organisms 
change their stability properties at will are yet to be established. 

6 Discussion and Conclusion 

We studied the locomotion, efficiency and stability of periodic swimming offish using a simple 
planar model of elliptic swimmer undergoing prescribed sinusoidal heaving and pitching 
motion in potential flow. We obtained expressions for the locomotion velocity for both small 
and finite flapping amplitudes, and showed how trajectories depend on key parameters, 
namely, aspect ratio 7, amplitudes A y and Ag and phases (fi y and (fig. Efficiency is defined 
as the inverse of cost of locomotion e. We studied the dependence of e on the parameters 
for both small and finite amplitude flappings. We observed that the efficiency maximizing 
parameters are approximately 3 < 7 < 4, A y w 1.3, Ag ~ ff,0y ~ (m + |)7r and (fig ~ nir, 
where n,m — 0, ±1,±2, ..., whose values are in excellent agreement with results based 
on experimental and computational motions of flapping fish, see [31 [161 H2] and references 
therein. 

We then studied the stability of periodic locomotion using Floquet theory. To our best 
knowledge, besides the work of Weihs which uses approximate arguments, this is the first 
work that rigorously studies the stability of periodic locomotion albeit in a simplified model. 
We focused on evaluating stability on the whole ((f> y , (fig) parameter space, and examined the 
effect of varying 7, A y and Ag on the resulting stability. We observed that stable and unstable 
regions in the (<fi y , (fig) plane evolved as these parameters change. Particularly noteworthy is 
the back and forth switching between stability and instability around the spots ((m+|)7r, nir) 
and (rwr, (m+ ^)tt). This switching is reminiscent to the observation in [10] that the motion 
of a heaving and pitching foil switches from coherence to incoherence and back to coherence 
when varying the aspect ratio of foil. In our study, we found a similar behavior when varying 
not only the aspect ratio but also the flapping parameters. This indicates that the swimmer 
can change its stability character by changing its flapping motion, and thus can easily switch 
from stable periodic swimming to an unstable, yet more maneuverable, state. Based on this, 
one could conjecture that, when it comes to live organisms, maneuverability and stability 
need not be thought of as disjoint properties, rather the organism may manipulate its motion 
in favor of one or the other depending on the task at hand. Clearly, this statement is 
speculative until verified by experimental evidence. To date, little is known experimentally 
on the stability of underwater periodic motions, let alone the stability of biological swimmers. 

Future extensions of the current work may include studying the effects of body deforma- 
tion and body elasticity on the observed stability of periodic swimming. In a planer potential 
flow, a model of articulated body that can undergo shape change, and with torsional springs 
attached to the joint to mimic the effect of fish muscle elasticity can be implemented (similar 
to the ones studied in [3 El E])- One may also be interested to study the effect of different 
frequencies of heaving and pitching motions (assumed to be equal in this work), and further- 

14 



more the effect of frequency in shape change (see [T7]) on the stability result. The current 
model makes the assumption of zero vorticity in the potential flow, which simplifies the anal- 
ysis and focuses on the added mass effect from the surrounding fluid. A natural extension is 
to include vortex shedding in the system, and use vortex sheet method [18] or Brown-Michael 
model [19] to study the effect of coupled body- vortex system on the aforementioned stability 
results. 



Appendix A: Hydro dynamic Forces and Moments 
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Figure 15: Exterior region of the ellipse of semi-axes a and b in the complex z-plane is mapped to the 
exterior region of a circle with radius r = (a + b)/2 in the £-plane. The mapping is given by (22). 



In potential flow, the fluid forces F x , F y and moment r can be obtained from the 
added-mass theory [20] or from the extended Blasius theorem [2Q1 [2TJ, [22]. In this appendix, 
we present both derivations and show their equivalence. 

The exterior region of the ellipse in the complex z-plane (z = x + iy) is mapped to 



the exterior region of a circle with radius r = (a + b)/2 in the £-plane, see Figure 15 The 
mapping is given by 



z - z 



i + j ) e l \ or equivalents, f = [Z !° )e - + ^ 



z — z. 



)2 e -2*0 _ 4c 2 



(22) 



where c = \/a 2 — b 2 /2. The complex potential of the fluid in £-plane is given by [22], 

Ur 2 - Uc 2 i8r 2 c 2 



W(t) 



s c 



e 



(23) 



where U = —z c e lB is the velocity of the mass center mapped into the £-plane. Therefore, 
the forces and moment exerted by the surrounding fluid on a moving body are given by the 
extended Blasius theorem [21]. In z-plane, 



F X T "Fy 



dwy , d 

— : — dz + ip— 
2 J dB \ dz J at 



ip 



i)B 



[z- z )——dz 
dz 



+ pA B z 



-Re 



2z (p [z — z j— — dz 
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dz 
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~dz~ 



dz + 
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- dw '. d2 



m 



dz 



(24) 
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where A B = ■nab is the area of the ellipse, dB is the boundary of the body, and the reader 
is reminded that the densities of the body and fluid are both p. Notice that the last term 
in r needs to be treated separately. All other integrals are analytic and, using residual 
theory, can be taken around an infinitely large circle instead of the boundary of the body, 
which greatly simplifies the calculations. For the last term in moment, since \z — z \ 2 is not 
analytic, one cannot use this technique. Instead, it needs to be integrated on the boundary. 



Substituting (22) and (23) into (24), one obtains the hydrodynamic forces and moment given 
by 



A , c 4 
F x = irp[ -'- — — + 2c 2 cos 29 ) x + 2npc 2 ysm29 - Anpc 2 9(xsm29 - ycos29), 






F y = np(- r + 2 ° - 2(? cos 29 j y + 2npc 2 x sin29 + 4npc 2 9(x cos29 + y sm29), ( 25 ) 
r = 2ixpc 2 (x 2 sin 29 - y 2 sin 29 - 2xy cos 29 - (?9 



Since 



2 m 2 -m 1 r 4 + c 4 m 1 + m 2 _, 4 
2vrpc = , irp ^ = , 2npc = J, 



the hydrodynamical forcing terms are equivalent with the expressions given in Q, which is 
repeated here for completeness, 

F x = - [—(mi + m 2 ) + (m 2 — mi) COS26 1 ] x -\ — (m 2 — m\)ysm.29 — (m 2 — mi) (x sin 29 — ycos29)9, 
F y = - [—(mi + m 2 ) — (m 2 — mi) cos 29] y -\ — (m 2 — mi)xsin29 + (m 2 — mi)(xcos29 + ysin29)9, 

r = —J9 H — (m 2 — mi) {x 2 sin 2$ — y 2 sin26 l — 2x?/cos2#) . 

And the governing equations are again repeated here 

m h x = F X} m b y = F y + F^, J b 9 = r + r^ p . (26) 

When expressed in a body-fixed frame, the hydrodynamic forces and moment take a simpler 
form in terms of the added mass coefficients mi,m 2 and J. Roughly speaking, as a body 
moves through potential flow, the body-fluid system behaves as an augmented body with 
modified mass and inertia that account for the added mass and added inertia due to the 
presence of the fluid. The added mass and inertia depend only on the geometry of the body 
and direction of motion. The Kirchhoff 's equations of motion in terms of the body-fixed 
frame variables are given by 

(m b + mi)Vi = —(mb + m, 2 )V 2 fi + F 1 , 

(m b + m 2 )V 2 = (m 6 + mi)Fifi + F 2 , (27) 

(J b + J)Cl = (mi - m 2 )V 1 V 2 + r fla P, 

where the body frame velocities and forces are given by 

Vx\ f cos 9 sin9\ fx\ a /FA / cos 9 sin9\ f F x 

V 2 J ~ \-sin9 cos 9 J \y) ' ' aM \F 2 ) ~ \-sin9 cos 9 J \F y 
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and r is the same form in inertial frame. Transforming (27) via a rotation 6 to inertial frame, 
one obtains the equations given in pi) and 0), and the hydrodynamical forcing terms are 
given by (El). It is then straightforward to verify that (26) and (27) are equivalent. 



Appendix B: Floquet theory 

For completeness, we repeat (fl8|) here 



where 



and 
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Hence, the equations of motion can be rewritten by taking inverse of 

q = g = M- 1 (f + F flap ). 
One can linearize above equation and obtain 

<9e 



/ ° \ 

P&ap 



\ r Aap / 



5q = J(£)<5q, where J(t) 

and the entries of the Jacobian JJ(£) are given by 
/ 2r 2 c 2 6sm26 2c 2 6(c 2 + r 2 cos26) 



dq 
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_ r 4 + c 4 
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-fi(3 



JJ13 
J 23 





(28) 



(29) 



(30) 



\ 



JJ14 

^24 
1 
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(31) 
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where 

a = xsm29 -ycos26, (3 = ysm26 + xcos29, fj, = 8r 4 c 2 /(r 8 - c 8 + 4r 4 c 4 ), 

2 r -i 2r 4 

F flap cos 29 - Ar 2 P Trf36 , J 14 = 4 _ 4 [-xr 2 sin 20 + y(c 2 + r 2 cos 20)] 



J 13 



J 23 



p7r(r 4 — c 



-n 



r* — c 1 



,2 



p7r(r 4 - c 4 ) 



•l 2c 4 
F flap sin20-4r 2 p7ra0 , J 24 = -^ j [x(c 2 - r 2 cos 20) - yr 2 sin2 



r -± _ c ^ 



The linear equation has fundamental solution matrix $(£), and any solution of (19) can 



be expressed as a linear combination of the four columns in $(£). Therefore, $(T) can be 



calculated numerically by integrating (19) with initial condition $(0) (it is convenient then 
to set $(0) = Z4X4), thus B = $(0) _1 $(T) can be solved numerically. Since x, F flap and 
r flap can k e so ] vec [ ^ mac hine accuracy, and analytical expression for J is available, B, and 
subsequently the characteristic multipliers, can be numerically solved to machine accuracy. 
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